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We present a new technique for time-domain numerical evolution of the scalar field generated 
by a pointlike scalar charge orbiting a black hole. Time-domain evolution offers an efficient way 
for calculating black hole perturbations, especially as input for computations of the local self force 
acting on orbiting particles. In Kerr geometry, the field equations are not fully separable in the 
time domain, and one has to tackle them in 2-1-1 dimensions (two spatial dimensions and time; the 
azimuthal dependence is still separable). A technical difficulty arises when the source of the field 
is a pointlike particle, as the 2-|-l-dimensional perturbation is then singular: Each of the azimuthal 
modes diverges logarithmically at the particle. To deal with this problem we split the numerical 
domain into two regions: Inside a thin worldtube surrounding the particle's worldline we solve for a 
regularized variable, obtained from the full field by subtracting out a suitable "puncture" function, 
given analytically. Outside this worldtube we solve for the full, original field. The value of the 
bkO' evolution variable is adjusted across the boundary of the worldtube. In this work we demonstrate 

^ ' the applicability of this method in the example of circular orbits around a Schwarzschild black 

, hole (refraining from exploiting the spherical symmetry of the background, and working in 2-f 1 

dimensions) . 

(N ; 

I. INTRODUCTION 

cj : 

I ' The motivation for this work stems from the problem of calculating the gravitational self-force (SF) acting on mass 
^ . particles in orbit around black holes. This problem has drawn much attention recently, in relation with the effort to 
model the inspiral of compact objects into massive black holes in galactic nuclei — one of the prime targets for LISA, 
the planned space-based gravitational-wave detector. The scientific merit from detecting such inspirals is potentially 
(N ■ high [H, but full exploitation of the gravitational-wave signal will require precise knowledge of the theoretical phase 
^ ] evolution of the waves, as predicted by general relativity, over a few years of inspiral. This, in turn, will require 
knowledge of the orbital evolution over a similar timescale, which, for sources of interest for LISA, would mean that 
the effects of the gravitational self force will have to be accounted for in the model. Thanks to the small mass ratio 
. characteristic of the relevant binary systems {fJ./M = 10~^-10~'', where jj, is the mass of the compact object and M 
is the mass of the massive black hole), the problem can be studied within the realm of perturbation theory, i.e., by 
considering the small perturbation caused by the SF to the orbit of a test particle moving in the fixed geometry of 
the central black hole. 
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[ There is now a well established theoretical framework for SF calculations in curved spacetimes 0, S 0, S @]j 
■ along with a practical calculation scheme for particle orbits around Kerr black holes 0, BQ- This method, dubbed 
"mode-sum scheme" , requires as input the local metric perturbation (or its multipole modes) near the particle, in 
a particular gauge — the Lorenz gauge. In the special case of a Schwarzschild spacetime (i.e., a non-rotating central 
5h ' hole) , the Lorenz-gauge metric perturbation can be obtained by solving the linearized Einstein equations directly, for 
each tensor-harmonic mode of the perturbation (using, e.g., numerical evolution in the time domain, as in p^). This 
method was applied recentl y fo r circular orbits in Schwarzschild, allowing a first calculation of the local gravitational 
SF for an orbiting particle |T]| . [Earlier calculations considered radial infall trajectories fl2| and static (supported) 
particles jl^, neither scenarios likely to be of relevance to LISA.] 

The work presented here is a first step towards tackling the problem in the Kerr spacetime. To obtain the Lorenz- 
gauge perturbation in Kerr, one may follow one of two possible avenues of approach. In the first approach, one 
first applies the Teukolsky formalism to solve for the perturbation in the Weyl scalars (fully decoupled into Fourier- 
harmonic modes if one opts to work in the frequency domain), and then reconstructs the corresponding metric 
perturbation in the Lorenz gauge. A procedure for Lorenz-gauge metric reconstruction is yet to be devise d. 33] . The 
alternative approach, which we pursue here, is to solve directly for the metric perturbation, as in Refs. [10lllll|. using 
time-domain numerical evolution of the Lorenz-gauge perturbation equations. This method offers a few important 
advantages: First, the problem of reconstructing the metric perturbation is avoided. Second, the behavior of the 
Lorenz-gauge perturbation near the particle, unlike that of the Weyl scalars, reflects the physical, isotropic form 
of the particle singularity, and is therefore more tractable. Third, the time-domain treatment best exploits the 
hyperbolic nature of the Lorenz-gauge perturbation equations. Finally, the time-domain treatment makes it easier to 
tackle particle orbits of arbitrary eccentricity. 

The main challenge in applying the above approach relates to the fact that the perturbation equations in Kerr 
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spacetime are not fully separable in the time domain. One can at most separate out the azimuthal dependence, then 
consider the evolution problem for each of the resulting 'm-modes', each of which being a field of 2+1 dimensions 
(2+1-D), depending on time and on 2 spatial coordinates. The non-separability of the field equations, on its own, 
does not pose a serious problem: Evolution codes for vacuum perturbations in 2+1-D have been developed and used 
successfully since the mi d-1990s [H, [H, [H [H [il, [il. The difficulty, rather, arises from the inclusion of a point 
particle as a source for the perturbation. Each m-mode of the resulting perturbation then diverges (logarithmically) 
at the particle, and accommodating this physical singularity on the discrete numerical grid becomes a major concern 
and the main challenge. 

Our goal here is to develop a scheme for handling the particle singularity within a 2+1-D evolution code. The 
idea is simple, and can be described as follows. As in Ref. [ll| (and unlike in [HI), we model the orbiting particle 
with a spatial delta-function distribution. The asymptotic form of the local perturbation field near the particle is 
then known analytically (e.g., [l^l). We then construct a function (given analytically) which (i) has the same local 
asymptotic form and (ii) is easily decomposed into m-modes. The difference between the full perturbation and this 
"singular" function defines a new, "regularized" field. The singular function is so designed that each of the m-modes 
of the regularized field is continuous at the location of the particle. We then use the (m-modes of the) regularized 
field as our numerical evolution variables. The full solution is simply the sum of the numerically-calculated regular 
field, and the analytically-given singular function. (Our "regular" function is, by construction, continuous, but not 
necessarily smooth; The regular function to be constructed in this work will have discontinuous derivatives at the 
particle's location. It should be stressed, in this regard, that our "regular" and "singular" variables do not necessarily 
correspond to Detweiler and Whiting's 'R' and 'S' fields the latter so defined that the 'R' field is a homogeneous, 
smooth solution of the perturbation equations.) 

To make it easier to control the global properties of our numerical evolution variable (especially its behavior at 
infinity and along the horizon), our evolution code will apply the above procedure only at the vicinity of the particle; 
far away from the particle it will utilize as an integration variable the full, original homogeneous field. To make this 
work in practice, we will introduce a reference "worldtube" around the worldline (in the 2+1-D space of the m-mode 
fields), whose "width" will be taken to be of order the background's radius of curvature, but will otherwise remain 
a control parameter in our numerical code. At each time step of the numerical evolution, the code will solve for the 
regular field inside the wordtube and for the original full field outside it, simply adjusting the value of the numerical 
variable across the boundary of the worldtube (using the known difference between the full and regular fields, being 
just the value of the singular function). In validating the numerical code, it will be important to monitor the amount 
by which the numerical solutions depend on the worldtube dimensions, and verify that this dependence becomes 
negligible with increasing grid resolution. 

The idea of representing a singular part of the solution analytically and solving numerically for the remaining 
regular part is reminiscent of the "puncture" method, often used in Numerical Relativity in representing initial data 



for spacetimes containing black holes [21[ . We shall call our singular variable a "puncture function" , and refer to our 



scheme as the "puncture method" , but we remind that here the idea of a puncture is applied in a different physical 
context. 

In this manuscript we demonstrate the applicability of the above method using a simple scalar-field toy model. 
To simplify the analysis still, we will consider circular orbits around a Schwarzschild, rather than Kerr, black hole. 
However, we will refrain from exploiting the spherical symmetry of the background geometry, pretending that the 
field equations cannot be separated into spherical-harmonic modes, and working in 2+1-D. The code we develop 
here should be expandable in a rather direct way to the Kerr spacetime and to eccentric/inclined orbits. We envisage 
applying a similar procedure for gravitational perturbations, but this would require much more development, including 
the formulation of the Lorenz-gauge perturbation equations is a format suitable for numerical evolution in 2+1-D. 

The structure of this paper is as follows. In Sec. [Til we decompose the scalar field in Schwarzschild spacetime into 
m-modes, and analyze the behavior of the individual modes near the particle, showing the logarithmic divergence. In 
Sec, mil we formulate our puncture scheme, select a particular puncture function, and analyze the asymptotic behavior 
of the regular field near the particle. In Sec. II VI we describe our 2+1-D numerical evolution code as applied for vacuum 
perturbations. We test it for numerical convergence, and check that, for initial perturbations with compact support, 
the late-time decay rate of individual multipole modes agree with that predicted by theory. We also test the 2+1-D 
vacuum solutions against solutions obtained with a 1+1-D code. In Sec.|V]we use the puncture scheme to incorporate 
a source term in our code, representing a scalar-charge particle moving in a circular geodesic orbit. We detail the 
numerical procedure in this case. Sec. IVII gives some results for the particle case, and presents a list of validation tests 
for the code. These include (i) test of point-wise numerical convergence, (ii) test of independence on the worldtube 
dimensions, and (iii) comparison with solutions obtained using a 1+1-D code. In Sec. IVIll we summarize this work, 
and discuss the application of our method for SF calculations. 

In passing, we briefly mention some related literature. Over the past decade, several authors have considered the 
evolution of black hole perturbations in 2+1-D, with or without a particle source. Krivan et at. [l3| wrote a 2+1-D code 
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to analyze the late-time power-law decay of homogeneous scalar field perturbations in Kerr spacetimes. Krivan et al. 
[l5| later examined also the late-time dynamics of the Weyl scalars associated with vacuum grav itational perturbations, 
by solving Teukolsky's master equation in 2+1-D. More recently, Pazos-Avalos and Loustofl^ presented an improved, 
fourth-order-convergent code in 2-l-l-D, for the evolution of vacuum perturbations of the Teukolsky equation. Particle 
orbits in Kerr have been tackled with a 2+1-D code by Lopez- Aleman et al. 16], Khanna [l3|, and, more recently, 
Burko and Khanna [2^ and Sundararajan et al. [19]. In these works (reporting on a series of improvements to the 
same 2-l-l-D code for evolution of the Teukolsky equation), the particle is represented by a smeared distribution of 
matter. The most recent of this works has achieved a reasonable accuracy in the far-field solutions, but the method is 
likely inadequate for accurate determination of the local field near the particle, which is essential for SF calculations. 
Sopuerta and Laguna [23| suggested the use of finite-element methods for an effective treatment of the particle. 
This idea (so far implemented for orbits in Schwarzschild |24]) shows much promise, but requires more development. 
Finally, Bishop et al. 25] have tackled the extreme-mass-ratio inspiral problem within the framework of full numerical 
relativity, i.e., by solving the full non-linear Einstein equations. This approach, too, requires more development. 
Throughout this paper we use metric signature = diag(— , +, +, +), and work in geometrised units, with G = c = 1. 



II. DECOMPOSITION OF THE SCALAR FIELD IN 2+1-D 



A. Physical setup: scalar particle in Schwarzschild 



We consider a pointlike test particle endowed with scalar charge q, moving in a circular orbit around a Schwarzschild 
black hole of mass M. In this work we ignore the SF, and assume the particle moves on a circular geodesic of the 
Schwarzschild background. Let Xp{T) denote the worldline of the particle (parametrized by proper time r), and 
introduce the tangent four- velocity = dXp /dr. Without loss of generality we work in a standard Schwarzschild 
coordinate system (t, r, 0, if) in which the orbit is confined to the equatorial plane, 6p = 7r/2. We then have 

< = [ipW-'^o -const, V2,t^ip(r)], ^ £/ fo[l, 0,0, lu], (1) 

where rp is the orbital 'radius', 

CO = d^p/dtp = (M/r3)i/2 (2) 

is the angular frequency (with respect to time t), and 

£ = ^upt = /o(l " 3M/ro)-i/2 (3) 

is the specific energy parameter, with /o = 1 — 2M/rQ. 

We take the scalar field $ of the particle to be minimally-coupled and massless. It then satisfies 

□<i>^-^(g"''$^V^)^^ = 5, (4) 

where g"'^ represents the background (Schwarzschild) metric, g is the background metric determinant, and the source 
term is given by 

J-00 V-9 

= =^^Sir ro)Sie - J)5(^ - ^tp). (5) 



r f 2 ' 



B. m-mode decomposition 



To reduce the problem to 2+1-D we decompose $ into azimuthal modes, in the form 



^ e'""^<i>"'{t,r,9). (6) 
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The individual m-modes are obtained through 

$™ = — / $e-'™'^dv'. (7) 

Note that, for future convenience, we take the principal values of the coordinate ip to lie in the range — tt < ip < t:. 
The scalar field equation ^ separates as 

= + + (/^ + 2r-i/'-) + 5"^ ($"9^9 + cot e - m^g-^-^^™ = S"" , (8) 

where the m-mode source reads 

gra ^ Z^ii _ 3M/ro)i/2j(r - ro)(5(0 - :^)e-™"*p. (9) 
Tq 2 

We note the relation ($'")* = (where an asterix denotes complex conjugation), which allows us to fold the 

TO < part of the sum in Eq. ^ over onto to > 0: 

00 

$ = + 2 ^ Re (e™-^*"*) . (10) 

m— 1 

To cast Eq. ([5]) in a form more suitable for numerical integration, we introduce the new variable 

$'" = r$'". (11) 

In terms of ^f"*, the field equation takes the form 

□m^m ^ _ J_ ^^m^ _^ cot - (2M/r + TO^ csc^ O) = -(/r/4)5", (12) 

where / = 1 — 2M/r, and m and ?; are the standard Eddington-Finkelstein null coordinates ('retarded' and 'advanced'- 
time coordinates, respectively), given by 

u = i — r*, ti = i + r*, (13) 

with 

r,=r + 2Mln('^l^]. (14) 



C. Behavior of '!>'" near the particle 

We now show that each of the modes diverges as a; — > Xp , and that this divergence is logarithmic (in the proper 
distance) . 

The singular behavior of the full scalar field near the particle is known to be described, at leading order, by (26l.[27| 

$(x) ~ ^. (15) 

Here x represents a point near the worldline, and e is the spatial geodesic distance from x to the worldline, i.e., the 
length of the small geodesic section connecting x to the worldline and normal to it. If x^ is a worldline point near 
X (not necessarily the intersection of the above normal geodesic with the worldline), and 5a;" = — Xp, then, at 
leading order in the coordinate distance, e is given by 

~ P^pdx'^Sxl^, (16) 

where Pap is a spatial projection operator reading 

PaP = QapiXp) + Ua{Xp)up{Xp). (17) 

Consider now a particular point x^ on the worldline, and let E be the spatial hypersurface t = tp, containing Xp. 
In the following we consider points x on S, and ask how $™(x) behaves Specializing to circular equatorial 
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orbits in Schwarzschild, we take, for simplicity (but with no foss of generality), ip-p = 0, and introduce local polar 
coordinates p, (f) in the i — plane: 



Sr 



= 9-90 = pr^^^ 



pCOS( 

psmc 



(18) 



1/2 

Then, at leading order, es = [j? + P^p^p(f'^^ , where the subscript S reminds us that e is evaluated on E, i.e, at 
t = tp. Substituting $ ~ q/es in Eq. we obtain 



$'"(p) 



27r 



-^dip (forx^Xp), 



(19) 



which describes the asymptotic behavior of <i>™ as one approaches the worldline along a t=const trajectory. Note that 
the particle limit corresponds to p ^ 0. 

To evaluate the above integral, we write it in the form 



. (p2 + p^^^2)l/2 



dip 



. (p2+p^^<^2)l/2 



dp 



1 



dip. 



(20) 



Using e 



l| = •\/2(l — cos mp) < m\(p\ (valid for \p\ < tt), the magnitude of the first integral on the right-hand 



side can be bounded, for any fixed value of p, as 



.(p2+P^^<^2)l/2 



dp> 



< 



m\p}\ 



. (p2 + P^^^2)l/2 



dp < 



1/2 



dp 



27rm 



P. 



1/2 



(21) 



(recalling P^^ > 0). Hence, this contribution is bounded at the limit p — > 0. Consider next the contribution from the 
second integral on the right-hand side of Eq. ([20]) . It gives 



dp) 



. (p2 + P^^<p2) 



1/2 



■In 



Po 



Pa 



(pI + p'V^'-po. 

0{p'), 



P 

2po 



(22) 



1/2 

where po = irPipip (depending on tq only), and where in the last step we expanded about /? — > 0. 
Collecting the results pT|) and ([22]), we conclude that, at leading order, 



*™(,^0) = -p-ln(^). 



(23) 



i.e., each of the m- modes of the scalar field diverges logarithmically with p, approaching the particle. (Note that p 
is the proper distance along geodesies in E emanating 'radially' from the particle.) Interestingly, the form of the 
leading-order divergence does not depend on the mode number m. 

Although we have restricted the above discussion to circular orbits in Schwarzschild, it is straightforward to repeat 
the analysis with an arbitrary point Xp along an arbitrary geodesic orbit in Kerr. The main conclusion holds in 
general: 'I'™ oc Inp as p — > 0, for any m. 



III. PUNCTURE SCHEME 



The divergence of $™ along the worldline is a serious concern when considering the numerical integration of the 
scalar field in 2+1-D. To deal with this difficulty, we introduce the following scheme. 

The scheme involves the introduction of a scalar field $p ('P' for puncture)^ given analytically, whose singular 
structure is similar to that of More precisely, we choose the field $p such that each azimuthal m-mode of the 
difference 



<I>R = $ - $p 



(24) 
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is bounded and continuous at the worldline. We then use the m-modes of $r, as new variables for the numerical 
integration in a region near the worldline. Specifically, we introduce a worldtube surrounding the worldline, the 
dimensions of which are kept controllable numerical parameters. Let T denote the interior of this worldtube, and dT 
its boundary. Let also <I>j^ and ^™ denote the m-modes of and $p, respectively: 

*R = 7^ $Re-""'^d^, $™ = -L ^pe-^'^d^. (25) 



27r7-^ 27r. 

The numerical scheme then utilizes the "regularized" variables $5^ for the part of the evolution which takes place 
inside T, while outside T it evolves the original fields <i>™. The value of the evolution variable is adjusted across dT 
using — $™ — $p . Thus, within our puncture scheme, the field equations to be evolved are 

= 5"" - 0"$^^' = S"™ in r, 
□m$m ^ Q outside T, (26) 

with = $™ - $™ on ar, 

where and $™ are given analytically. Of course, once the continuous fields $p are solved for, the full scalar-field 
modes can be simply constructed through = 4>J^-|-<i>j?. Note that, depending on the form of the puncture function 
$p , the source may have support anywhere inside T (not necessarily confined to the worldline) . 



A. Choice of the puncture function $p 

We wish to construct a function <i>p which (i) reproduces the singular behavior of the full field $ at a; — > Xp; (ii) is 
sufficiently regular away from the particle; and (iii) is easily decomposable, in analytic form, into m-modes. 
Consider the puncture function 

<^>p{x;xp) = ^, (27) 
ep 

with 

ep = ^Jp^ + 2P^^{1- cos Sip) . (28) 

Here, as in SecHH Ja;" = a;" — Xp, Pap are tensorial coefficients as defined in Eq. (fTT)) . and p [same as in Eq. p8|) ] is 
given explicitly by 

p^ = P„Sr^ + PggSe^ . (29) 

Note that here we do not regard the coordinate differences Sx as necessarily small. The definitions in Eqs. (|27p - (P51) 
are taken as exact, for any value of 6x within the worldtube T. 

Since 2(1 — cosScp) — Scp'^ + 0{6(p'^), the function ep coincides, at leading order in 6x, with the function e [see Eq. 
(fTB)) ]. evaluated on the hypersurface t = tp. Therefore, at leading order in Sx, the puncture function <l>p coincides 
with ^{t = tp). As desired, $p is singular only at the location of the particle {p = Sip — 0), and is regular (C°°) 
anywhere else. Finally, as we show below, our function $p is easily decomposed, in explicit form, into m-modes. [s^] 



B. Continuity of the modes $r 

We now show that, with the puncture function selected above, the modes = $™ — are finite and continuous 
for all r and 6. Since both $ and $p are regular {C°°) away from the particle, then so is $p, and so are its modes 
We hence focus on the behavior of the modes $p at the location of the particle {Sr = 56 ^ 0, or, equivalently, 
p = 0), aiming to show that they are C° there. 

For this discussion, we will need to consider higher-order terms in the asymptotic formula (jlSp . ft was shown in 
Ref. [i^l (by considering the Hadamard expansion of the retarded Green's function for the scalar field) that, near the 
particle. 



^x) = ^ + h{x), 



(30) 
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where /i is a C" function (i.e., continuous but not necessarily differentiable) . The spatial geodesic distance e can be 
expanded in terms of the coordinate difference Sx"' in the form 

^el + Qap-yixp)Sx°'5x^Sx'' + 0{Sx^), (31) 

where Eq = Pai3{xp)Sx°'6x^ [the leading-order form, as in Eq. (|16p]. and Qafij are certain coefficients depending only 
on Xp (they are given exphcitly, e.g., in Ref. ^26]). Substituting for e from Eq. (I3T|) . Eq. ([30]) becomes 



^ - - ^gQo/a /"^ ""3 ""^ + /2(a:), (32) 
where /2 is C°. 

In what follows we fix Xp, and consider the behavior of $ (and <l>p) on the hypersurface t — tp, denoted E as before. 
Recalling the definition of ep in Eq. we have, on S, 

el = el + 0{5x^) (33) 

and hence (near the particle) 



1 _ 5x'^5xf^5x'< 



Thus, on S, 



$P = — + 0(fe). (34) 
eo 

1 5x°'5x^5x'^ 

$R = $-$P==--gQa/37 3 + /3(a;), (35) 

where /a is yet another C*^ function. In this expression Ja;" — dr, SO, or and 

eo-(p'+P^^V)'/', (36) 

where, recall, p is given in Eq. (|^^ . We now write <I>j^e~*'"'^''' = 'I'Ril + 0((5(^)] (for any m, at small l^i^sj), and notice 
that, by virtue of Eq. ([35|) . the contribution to ^pe"*""'*' from the 0((593) terms vanishes as eo — > 0. Hence, this 
contribution is C", and we may write 

■ I 1 5x°'Sx/^Sx'' 

^^^-^m5^ = -^qQ^,, 3 + /r(^), (37) 

where /"(a;) is a C° function (depending on m), and where the first term on the RHS is m-independent. For simplicity 
(but without loss of generality) we take ipp = 0, giving Sip — ip. The m modes of ^p, are then given by 

'^'l = -\qQ.pJ V d^ + fE^, (38) 



where the integral f^{r,9) = J^^f"'dip is necessarily a C° function of dr and 69 (since the integrand /™ is a C° 
function of dr, 59 and Lp). 

It remains to show that the term c>c Qap-^ in Eq. ([55]) is . (Note that the integrand in this term is not necessarily 
C'^ .) To this end, we write {~q/2)Qai3'ySx°'Sx^Sx'' explicitly as a polynomial in ip, in the form + P2'P + Pi'f'^ +Pof^ j 
were Pn{Sr, SO) are polynomials in Sr and 69, each of the form X]fe=o 0'k6r^ 69"^^^ (with at constant coefficients). The 
contributions from the terms cx P2,Pa to the integral in Eq. (j38p vanish from symmetry The contribution from the 
term cx ps reads 

(X Am r 27rp3{6r,60) 
P3{5r,60) . „ „ 5-^ = — — ^===^0 39 

as p — > 0. The contribution from the term cx pi also vanishes at the limit p 0: 
Pi{6r,60) 



V^dif Xa\V-3/2 1 P0+ VPo+P'^\ ( 2 , 2N-1/2 

.(TTTV^ = 2p,(fr,5^)F^/ ^ln(^ J-po(Po+H 

Hence, the integral in Eq. ([55]) vanishes as p ^ and is therefore a C° function of 6r and (50. 

The above verifies that the modes are each continuous at the location of the particle (and elsewhere). We 
do not expect, however, the derivatives of to be continuous. [That the derivative are likely to be discontinuous 
is suggested, for example, by the form of the contribution evaluated in Eq. (j40p .] In the numerical scheme to be 
developed in Sec.|V]we shall assume explicitly that the solutions $p are continuous. 
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C. Expressions for the puncture modes $5? 

To implement the puncture scheme set out above [Eq. (j26[) ]. we need exphcit expressions for the puncture modes 
<i>™ and for the regularized source modes S^. We start by obtaining the necessary expressions for <i>p . 
With the above choice of $p [Eqs. (P7|) - (P9)l ]. the modes <i>p are given by 



27r 



-TT y/p'^ + 2Pyy(l - cos Sip) 

dx 



J-TT-ujt^ \J + 2Pipip(l - cosx) 

<? ^-^,y^u^U T COs(mx) 



27r ^p2 + 2P^^(l-cosx) 

where in the second integral we have changed the integration variable QSip^x = 8Lp = ip — wtp, and where in the 
third integral we (i) shifted both integration limits by Lot-^ (noticing the integrand is periodic with period 27r), and 
(ii) made use of the fact that the imaginary part vanishes by symmetry. For all m = 0, 1, 2, . . ., the last integral can 
be represented in terms of complete elliptic integrals. We find 

„ —imuitp 

*P = TJT- bS(p)7^(7) +P^(p)7i?(7)] , (42) 

where 

7^[l + pV(4P^^r'/', (43) 

K{j) and E{j) are complete elliptic integrals of the first and second kinds, respectively (as defined in [1^), and 
and are certain polynomials in p^. We tabulate these polynomials in Appendix \X\ for m = 0-5. 

D. Expressions for the source modes 

Within our puncture scheme, the source for the field $r inside T is S'r = S — and its to- modes are given by 

^o) = ^fiS- D^p) e-""^d¥.. (44) 
With the function $p defined above, and using (p = 6(p + cjtp, this takes the form 

= :f e-™'"*^ (5i/r + S2IT + Ssir + SdT) , (45) 

ZTT 

where the /™ (depending on rg only) are the integrals 

/{" = / ep^/^e-™'^'^d(^^), 

J — TT 



— TT 



= / ep^/^sin2 5(^e^""'''^d((5v3), (46) 



and where the Sn are m- independent functions of r and 9 (as well as ro), given by 

51 = Prrf{r) + 2r-'^Prr{r- M)5r + r-'^Pgg{l + 5ecoii 

52 = P^^ [r-^ sin-2 e - coV/(r)] , 

53 = -3Plfir)6r^^3r-^P^g5e\ 

,2 



54 = -3P^^ [r-2 sin-2 6 - c.V/(r)] . (47) 
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In obtaining Eq. (|45p from Eq. (|44)) one should note the following: Firstly, the function $p depends on t, through 
Sip = tp — ipp = ip — Lutp = tp — ujt (as in our construction we take t = tp)] this should be be taken into account properly 
when evaluating Secondly, the source Sb, contains no distributional component (i.e., the delta functions in S 

and □'fp "cancel each other" precisely) [HI . 

The integrals /{" 4 can once again be expressed in terms of complete elliptic integrals. Introducing the dimensionless 
"local distance in the r-9 plane" , 

(48) 



jl/2 

[in terms of which the quantity 7 of Eq. (I43p reads 7 = (1 + /5^)^^/^], we have 

C=i,2 = ^^"^'/'7fc(p)i^(7) + r'Cij(/5)i?(7)], 

IT = P^^'^l WAmil) + ~P~hnE{~P)E{l)] , (49) 

where p™^^ and p™^ are polynomials in , which are all non-zero at p — > 0. (Note the particle limit corresponds to 
p ~f 0+, or 7 ^ 1^.] The polynomials p™^ and are tabulated in Appendix VK\ for m — 0-5. 



E. Behavior of 5*^^ near the worldHne 



Even though the modes are finite (and continuous) at the worldline, the source modes 5^ can still diverge there. 
Indeed, as we show below, diverges like as p ^ (with coefRcient that depends on the direction of approach in 
the r-9 plain). This is a serious concern when it comes to numerical implementation, since the finite-difference scheme 
would normally require to evaluate the source S'J^ also at p = 0, where it diverges. We will deal with this complication 
by integrating "by hand" numerical grid points which lie on the worldline (this procedure will be described in Sec. 
|V|) . For this, we shall need some information on the asymptotic form of near p — Q. It will prove necessary to 
have at hand the form of 5^* up to 0{p'~') (inclusive). We now derive the necessary asymptotic formula. 

Consider the form of as given in Eq. (^5)1 . The functions Sn are easily expanded in powers of Sr and S9. We 



wish to rewrite this expansion in terms of local polar coordinates as in Eq. psp . However, our numerical coordinates 
will be based on t, rather than t, r, and it will prove advantageous to replace p as our local polar coordinate with 
a new coordinate, based on 5r^ = r» — r^Q. We hence introduce the new local polar coordinates p*, 0*, defined by 



(5r* = p»/o ^ cose/)*, 
S9 — p^TQ^^sin^,. 



(50) 



(The coordinates p*,0* coincide with p, at leading order in p, but deviate at higher order.) Using 6r — fodr^ + 
^/'(''o)'^'''* + • • • and Eq. ([50]) we can then express each of the 5„ as an expansion in p,. To expand the /™ in powers 
of p*, we first expand the Elliptic functions in Eq. (I49|) in powers of p, using 



X(7) 



1 



ln(x/4) - 7X'[ln(x/4) + 1] +0(xMnx), 



^^(7) = l--x'[21n(x/4) + l]+0(x'lnx), 



(51) 



[Eqs. (8.113-3) and (8.114-3) of Ref. ^] where x = ^1 - 7^ = p(l + p^)"^/^ (recall the particle hmit corresponds 
to p ^ 0+, or 7 — > 1^, and so also to x ^ 0+). We then re-expand /™(p) in powers of p, using 



M 



2 £3/2 



M(llM-4ro) 



: sm 



2M 
3^ 



pI + 0{pt), 



(52) 



obtained by substituting 6r — /o<5r* -I- ^f'{ro)5rl in Eq. ([2^ . and then substituting for Sr^, and S9 from Eq. ((50)1 . 
Inserting all the above expansions in Eq. (j45p . we obtain the following asymptotic formula for the source modes: 

aicj). 



Qm ^ —imujtp 



+ Anln(P*/4)+/3" 



+ ©(p* Inp*), 



(53) 
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1/2 

where = p^^K^P^pip ), and the coefRcients are given by 



8(1 - Af/rp) 



cos (/), sin^ (f). 



(54) 



An 



1 



(55) 



/3"(0*) 



/3™ + ;9i cos 24>^ + ^2 cos ^^ + /^a cos 60, . 




(56) 



In the last expression /3™ are constant coefficients (depending on tq and m only), of which only will be needed in 
what follows — this coefficient is given explicitly in Appendix |B] for m = 0-5 . 

We note the following: (i) At leading order we have oc with coefficient that depends on m only through the 
trivial factor e"'™"^*?. (ii) The leading-order divergence of is direction-dependent (it depends on the azimuthal 
angle in the r-9 plane), (iii) Upon averaging over directions, the ^ p^^ divergence of cancels out; the direction- 
averaged singularity of S'^ is only (x Inp, . We will make good use of this latter observation below, when setting out 
our numerical scheme. 



In this section we develop our 2-f 1-D numerical evolution code and test it for vacuum perturbations. To this end 
we set, for now, S"'" = in the field equation (fT2|) . and consider the vacuum evolution of prescribed initial data. For 
each mode to > we discretize the field equation using a 2nd-order-convergent finite-difference scheme, on a fixed 
2-1- 1-D mesh which is based on mixed characteristic and spatial coordinates. We test the validity of the code by (i) 
demonstrating 2nd-order convergence, (ii) examining the late-time decay pattern of compact initial perturbations, 
and (iii) comparing the numerical solutions to those obtained using evolution in 1+1-D. 



FIG. 1: A diagram illustrating the 2-l-l-D numerical domain. The grid is based on characteristic (Eddington-Finkelstein) 
coordinates u and v and Schwarzschild coordinate 6. Initial data are specified on the (null) surfaces t; = and ti = 0. Boundary 
conditions are specified at the "poles", 6 = Q,n. 

The 2-l-l-D numerical domain consists of a "stack" of staggered double-null 1+1-D grids, each based on w, u 
coordinates — see Fig. [TJ We denote the grid spacing in each of u,v by h, and the grid spacing in 9 by A. The 
evolution starts with initial data on the two hypersurfaces defined hy v — vo and u — uq. (In the circular-orbit case 
considered in the next section, these will be taken such that the initial vertex uq,vo corresponds to r = rg, t = 0, 
where r = tq is the orbital radius.) The numerical evolution proceeds first along 9, then along u and finally along 
V. That is, for each v value we solve for all u, and for each it, v values we solve for all 9. Boundary conditions (see 
below) are placed along the two surfaces 6* = 0, tt, representing the two polar axes. 

As the grid is not based on purely characteristic coordinates, we must constrain the relation between h and A. On 
theoretical grounds, for the scheme to be stable it is necessary that the numerical domain of dependence contains the 
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A. 



Numerical domain 



hx h 
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physical, continuum domain of dependence at each point in the evohition ("Courant condition"; see, e.g., [231). In our 
scheme, a grid point at (t, 6) will effectively require data from points (< — ft,, ± A), so the above condition translates 
to A/ft > f^/'^/r. The function f^^'^/r attains a maximum value of ~ 0.19245 il/^^ (at r = 3Af), so to make sure 
that the Courant condition is met everywhere, we shall always take A/ft > 0.2 M~^. 



Finite difference scheme 




FIG. 2: A single numerical grid cell, of coordinate dimensions h x h x (2A). Finite difference approximations are made about 
the point c in terms of points 1-8. The equations are then rearranged to give an evolution scheme for the field at point 1, based 
on the values at points 2-8 solved for in previous steps of the evolution. 



Our numerical evolution scheme is constructed from finite-difference approximations to the terms in Eq. (jl2p . 
centred on the point c as shown in Fig. [21 At each stage we solve for the point 1 based on information from the points 
2-8. We then rearrange the resulting formula to obtain an evolution scheme for the point 1. The various terms in Eq. 
(I12p are approximated at point c using the centred formulas 



ft2 



O(ft^), 



(57) 



*5 



2(5-™ + -qjf) 



2A2 



o(A^ft2), 



(58) 



4A 



(59) 



0(ft2). 



(60) 



All r, ^-dependent coefficients in the field equation are simply evaluated at point c. Solving for we obtain our 
finite difference scheme for the vacuum case: 



Sri 



2^-^ - 2^-^) /A^ 



- cot 6»e (4-^ + -^'^ - ^r™ - ^r™) /(2A) 

- (2M/re + csc^ 9^) (^-^ + ^'™)] + ©(ft^A^, ft-*). 



(61) 



Here Tc and 9^ are the values of r and 9 at point c, and /c = f{fc)- 

For a fixed ratio A/ft, the finite-difference scheme ([M|) has, effectively, a local discretization error of ©(ft"'), leading to 
a global (accumulated) error of 0fft^).[36| Hence, we expect the algorithm to exhibit quadratic point- wise convergence 
(for smooth initial data). 
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C. Boundary conditions at the poles 

The boundaries in u, v are null and thus are never encountered during the evolution. On the other hand, at the poles 
(0 = 0,7r) we require suitable boundary conditions. One can obtain the necessary conditions by imposing regularity 
of the field at the poles: Each azimuthal mode m 7^ has harmonic dependence on ip; continuity of the field across 
the poles (where if is indefinite) therefore implies that the field must vanish there. The remaining, axially-symmetric, 
m = mode is symmetric across each pole (invariant under ip — > and so for the field to have continuous 

derivatives there, these derivatives must vanish. The physical boundary conditions at the poles are therefore 

^m#O(0 = 0, tt) ^ 0, *;r°(^ = 0, ^) = 0. (62) 

To implement these conditions in our code, we simply set 4* = at the poles for all to 7^ 0, whereas for m = we 
use the extrapolation 

^m=0(g, ^ 0) = i [4*"="(6' = A) - *"=0(6i = 2A)] + ©(A^^), 

^m=o^Q = tt) = i [4*"="(e' = TT - A) - *'"="(6i = vr - 2A)] + 0{A'^). (63) 
[Here the error term is 0{A'^), rather than 0{A'^), since v[/™=o is an even function of 9 at the poles.] 

D. Tests of vacuum code 

For the following tests we specified initial data in the form an "outgoing" narrow pulse starting at vq,uo, with 
a certain 6'-profile chosen differently for each of the tests (see below). In all cases we took vq = r*(r = 7M) 
and uq = —r-j,{r = 7M). We selected A such that, at the lowest resolution, tt/A is an integer number (and so 
an integer number of A intervals fit into our grid between the two poles). In all cases we fixed the ratio A/h 
at 27r/5M~^ ~ 1.26 A/^^. This is safely above the Courant limit, and for our lowest resolution {h = M/A) gave 
sufficient ^-resolution (A = tt/IO) to resolve the lowest few multipoles. 

1. Numerical convergence 

We tested the point-wise self-convergence rate of the above scheme by examining the solutions along various 1-D 
cross sections of the 2-l-l-D grid, for a geometrical sequence of decreasing h values approaching /i ^ 0. Specifically, 
we looked at *™(i) along (r, 9) = (JM, 7r/2) and at ^'"(6') along (t, r) = (400M, 7Af), for resolutions h = M/2, M/4, 
and Af/8. For initial data, we took a narrow distribution, centered at {v,u,9) = (t^o, uq, 7r/2). We deliberately chose 
a discontinuous initial distribution (a narrow square pulse), which simulates the situation in the particle case (see 
below) and allows us to assess the effect of non-smoothness in the initial data on the convergence rate. 

As demonstrated in Fig. [31 the vacuum numerical evolution shows a clear second-order point-wise convergence at 
late time. Early in the evolution, multiple reflection of the discontinuous data off the 6'-boundaries introduces large 
numerical error [for grid cells in which ^ is discontinuous, the finite-difference formula ([61]) has a local error ofO{h^) 
rather than 0{h'^)]. However, over time the discontinuity dissipates, and quadratic convergence is retained. 

2. Late time tails 

In theory, after the initial burst of radiation and ringing phase, the field should settle down to a power-law decay 
at late time. The exponent of this decay is determined predominantly by the lowest multipole number £ present in 
the data, but is otherwise independent of the shape of the initial data . In the case of a Schwarzschild background 
(where different i!-modes do not couple) and compact initial data with angular dependence of a pure £-mode, we 
expect late-time tails of the form ^I^™ cx t~'^^~^ at fixed r (i.e., for t ^ M with t ^ r), and of the form '5'" cx 
at null infinity (i.e., for v ^ M with v ^ u) i3C| . 

To test these predictions with our vacuum code, we specified initial data in the form of a compact outgoing 
pulse of a pure ^-mode content. Specifically, we took 5'™(u = uq) = for all v,9, and \I'"'(w — vq) — sin^[7r(u — 
uo)/{8M)]Pirn{cos9) for < m < 8A/, with ^/"^{v = vq) ^ for u > 8M. Here Pi„i{cos9) is the associated Legendre 
polynomial. The left panel in Fig. [H shows the late-time decay tails of our vacuum solutions at fixed r(= 7Af), for 
£ ~ TO = 0, 1, 2. The right panel in Fig.[51shows the decay tails as a function of retarded time u at large v{= lOOOAf), 
approximating null infinity. The decay rates are in excellent agreement with the theoretical prediction. 
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FIG. 3: Numerical convergence test for vacuum perturbations. Upper, middle and lower panels correspond to m = 0,1,2, 
respectively. Each of the 6 panels labeled 'm = . . .' displays jX the relative difference = |(*(^) - - *(8))|> 

where "Jj^) is the solution obtained with resolution h = [5Af/(27r)]A = M/n. A value of unity indicates quadratic convergence. 
Left panels show along {r,6) = (7M,7r/2); right panels show jS<S)^^i{e) at {t,r) = (400M,7M). For S'f^iie) we also 

show the solutions themselves (3 small panels): ^'j^j in dotted line, '^^^l) in dashed line, and "i!^) in solid line. These 

serve to demonstrate how the late-time solutions are dominated by the lowest allowed ^-mode for a given m, i.e., imin = \m\. 
The solutions show good 2nd-order numerical convergence at late time. At the early stage of the evolution, the solutions are 
affected by the initial-data discontinuity, bouncing back and forth between the two poles; however, this effect gradually dies off 
over time through dissipation. The noise in <5*™r^ at very late time is due to round-off truncation error, which kicks in when 
the amplitude of drops very low. 

3. Comparison with 1+1-D solutions 

The best quantitative test of our code comes from comparison with results obtained using an independent evo- 
lution code formulated in 1+1-D (timc+radius). For this comparison we wrote a 1-|~1-D code similar to the one 
developed in Ref. [3ll |. In the 1+1-D treatment we construct $ through an expansion in spherical harmonics, 
$ = (l/r)y^^^ ^^"'(t, r)y^™(6>, (jg), where the time-radial functions 5'^™(t,r) are obtained using characteristic evo- 
lution in 1+1-D (see [31| for details). Suppose that is a 1+1-D solution for given £,m and for initial data in 
the form of a compact outgoing pulse with some profile U{u). Suppose also that is a 2+1-D solution with the 
same m, for initial data in the form of a compact outgoing pulse with a profile U{u)Pim{cos9). Then, we expect 
the solutions to be related by \E'^]^(t, r, 0) — ai,n^^{!^i{t,r)Pim{cos6), where agm are the normalization coefficients 
appearing in the relation between the spherical harmonics and the Legendre polynomials: Y^"^ = agmPem{cos 9)e'^™''^ . 

For the comparison, we ran the 2+1 code with the same initial data as for the tail-test above. We then ran the 
1+1-D code with corresponding i,m and w-profile. The plots in Fig. [5] display results from this comparison (for 
£ = m = and i = m = 1, along lines of constant r,d). We find a good agreement between the 2+1-D and 1+1-D 
solutions. 



V. NUMERICAL IMPLEMENTATION: CIRCULAR ORBIT 

A. Inclusion of the particle; the worldtube 

We now come to the main part of our analysis: the inclusion of the particle through the puncture scheme described 
in Sec. IIIIl We consider the physical setup described in Sec. Ill Al i.e., a scalar-charge particle set in an equatorial 
circular geodesic orbit with radius r — vq around the Schwarzschild black hole. The scalar field equation now has 
source S, given in Eq. ([5]). In our 2 + 1-D numerical domain the particle traces a straight line along u = u + 2r,o, 
9 = tt/2, where r^,o = r,(ro) (see Fig. (5]). We set up the grid such that the initial vertex {vo,uo) corresponds to t = 
and = r,o; namely, we take vq — r,o and uq = — r^o- We select the 9 grid separation A such that tt/A is an even 
integer. With this setup, which turns out most convenient, the particle cuts straight through grid points precisely 
every At — h. 
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FIG. 4: Late-time power-law decay tails of with initial perturbation made of a single, pure £ = m mode. Left panel 
(solid lines): tails along lines of constant r = 7M,0 — 7r/2. Right panel (solid lines): tails at "null infinity", read off along 
V — const — lOOOAf, again at ^ = n/2. (For clarity, part of the ringing phase data has been removed in these figures.) The 
dashed lines are reference lines a t~^^~'^ (left panel) and oc u~^~^ (right panel), showing the theoretical asymptotic slopes. 
The gradual deviation from the predicted slopes in the "null infinity" data is explained by the fact that the large-u regime (m 
comparable to v) no longer approximates null infinity. 




FIG. 5: Comparison between vacuum solutions obtained with our 2+1-D code, and solutions obtained independently using 
1-l-l-D evolution. Both solutions correspond to the same physical initial data, containing a single, pure I, m mode (see text for 
details). We compare them here as functions of t, at fixed r = 7M and 6 — -k 12. The left and right panels display ^ = m = 
and £ — m = \, respectively. The upper panels show, superposed, both *I'5^i(t) (dotted line) and *I'S+i(t) (dashed line). The 
lower panels show the relative differences 2 |(^'i!^2 — ^r+i)/(^i+2 + ^i+i)!- The small relative difference is (presumably) due 
to finite-differentiation errors in both codes, which in the 2-l-l-D code also include a small amount of "contamination" from 
coupling to higher £-modes. (For ^ = m = the tiny relative difference is dominated by noisy round-off error.) The good 
agreement between the 2-l-l-D and 1-l-l-D solutions provides a strong vaUdation test for the 2-l-l-D code. 



To solve the sourced evolution problem, we implement our puncture scheme as formulated in Eq. (|26p . We first 
introduce a "worldtube" T within the numerical grid. For convenience, we choose a worldtube with a uniform 
rectangular cross section: For any fixed time t we take it to be the region ~ ^r, /2 < < r*o + <5r» /2, 7r/2 — 5^/2 < 
B < 7t/2 + 60/2, where the "width" 6r, and "height" 6$ of the tube are kept as (two independent) control parameters 
in our analysis; See Fig.[6]for an illustration of the worldtube setup. We will typically take Sr, and roSg to be of order 
a few M. Among the robustness tests for our code, we will establish that the numerical solutions are independent of 
6r.^ and 6g (up to numerical error which decreases with grid size). 
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FIG. 6; The worldtube configuration. The sketch illustrates the geometry of the numerical domain in the circular-orbit case. 
We show a portion of the numerical grid, containing the worldline (dashed line), and the worldtube surrounding it (shaded 
volume). More details are given in the text. 



B. Finite-difference scheme 

For our numerical treatment, we reformulate the puncture scheme (j26[) in terms of the field variable = r^™, as 
in the vacuum case. The scheme becomes 

= -(/r/4)5™ = in T, 
□m^m ^ Q outside T, (64) 

with = ^r™ - r$]? on dT , 

where, recall, the operator is defined in Eq. \12\ . and SJJ* and <i>™ are given analytically in Eqs. ([l5|) and (|42p . 
respectively. The evolution algorithm is similar to the one applied in the vacuum case: Starting with initial data 
ow V — vq and u — uq (see below), we integrate along "planes" of fixed w, where on each such plane we integrate 
along "lines" of fixed u. Consider again a typical, single grid cell as depicted in Fig. O The values at points 2-8 are 
assumed to have been solved for in previous steps, and we need to obtain the value at point 1. The algorithm first 
labels each of the points 1-8 as either 'out' or 'in', depending on whether it lies outside or inside T, respectively. If a 
point lies on dT it is labeled 'in'. Four cases are possible: Case 1: All points 1-8 are 'out'. In this case the integrator 
implements the vacuum scheme to solve for ^f™, just like in the global vacuum case. Case 2: All points 1-8 
are 'in'. The integrator then implements a different scheme, described below, which is based on the sourced equation 
'-'*^R — ■ Case 3: Point 1 is 'out', but some of the points 2-8 are 'in'. In this case the values of the 'in' points 
are adjusted according to ^™ = ^tj^ + r^j?, after which the integrator solve for 5"™ using the vacuum scheme 

(ICTj) . Case 4: Point 1 is 'in', but some of the points 2-8 are 'out'. Then the values of the 'out' points are adjusted 
according to vp™ '^'^ = ^™ — r<I>™, after which the integrator obtains at point 1 using the sourced-equation 
scheme described below. This way, the algorithm effectively solves for 5'™ outside T and for 'i'^ inside T, adjusting 
the integration variable across the boundary dT. 

We now describe the finite-difference scheme applied inside T. Referring again to Fig. [51 we consider the case where 
point 1 is labeled 'in', and assume the value of ^'r at points 2-8 has been obtained in previous steps (possibly through 
the adjustment ^i™ ^ ^f™ = ^f™ - r$j?). The field ^-jj obeys the inhomogeneous equation □^^'j:^ = Z^, where the 
source term is known analytically, and has a definite finite value everywhere, except on the worldline. To obtain 
^'r at point 1, we first write finite-difference approximations for □^^E'^^ centered at point 0, as in Eqs. ([57]) -([60 ]) . 
If point 1 is off the worldline, then so is point c, and we include the source term by just evaluating Z^ at point c. 
Solving for yields the finite-difference formula 

^-^^ = [RHS of Eq. Unil), with ^'J^ ^ + /i^Z™ (in T, off the worldhne), (65) 

where Z^^ is the value of Z^ at point c. This scheme has local discretization error of 0{h'^), with a global (accumu- 
lated) error of 0{h'^), just like the vacuum scheme. 

The source S^, and so also Z^ = — (/r/4)S'^, diverges at the worldline in a manner described by the asymptotic 
formula (|53p . i.e., like ^ p~^, with amplitude depending on the direction of approach in the r-0 plane. The divergence 
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of the source poses a technical problem when it comes to numerical implementation: Even if we re-arrange the grid 
such that grid points are always avoided by the worldline, still the rapid growth of the source near the worldline would 
be difficult to accommodate numerically. (The problem will show up more acutely for non-circular, non-equatorial 
orbits, where it will be more difficult to assure that the worldline does not "come too close" to any of the grid points.) 
A natural solution to this problem could be achieved within a higher-order puncture scheme, in which the puncture 
function is taken to account for additional, subdominant terms of the local field. We leave the formulation of such 
advanced scheme for future work; here we will continue to use our leading-order puncture, and demonstrate that even 
this simple scheme can yield numerically-robust solutions. 

Rather than trying to "avoid the worldline" with a suitable layout of grid points (a strategy which will not be useful 
anyway for more complicated orbits), we take here the worldline to cross straight through grid points. To derive the 
finite-difference scheme for points crossed by the particle, we integrate the field equation locally "by hand", as we 
describe in what follows. This guarantees that need never be evaluated at a distance smaller than Ar^, = h/2 
from the particle. We envisage applying a similar local-integration procedure for generic orbits, which would save the 
need to carefully lay out a "particle-avoiding" grid. 



C. Treatment of worldline points 

Referring, once again, to the grid cell illustrated in Fig. [21 we consider the case where point 1 (and so also points c 
and 4) lie on the worldline. We assume that the value of ^'J^ at points 2-8 has been obtained in previous steps, and 
we need to approximate the value at point 1. Inside the cell, the source term diverges as oc p~^. The integaral 
of Z^ over the volume of the cell should therefore yield a finite value. Moreover, we observe in Eq. ([55)) that the 
direction-dependence of the leading order, oc divergent term of Z^ is such that the contribution from this term 
vanishes upon integrating over all directions. This suggests a method for deriving a finite-difference formula for points 
on the worldline: Based on the values of at points 1-8 of the grid cell, write a finite-difference approximation for 
the integral equality 

/ {a^^^)dV^ f Z™dy, (66) 

Jccll Jcoll 

where the integral is evaluated over the 3-D volume of the grid cell shown in Fig. [21 and dV — du dv d9 is a coordinate 
(not proper) volume element; Then solve the resulting discrete algebraic equation for the value of at point 1. 
To formally discretize the LHS of Eq. ([66]) , we write 

/ (a^i^^) dV = 2A (M/^i - - *fl3 + *S4) 

J cell 

-/i2A(2Af/ro + m2)(*^'2 + *™3)] +0{h^\nh). (67) 

Here 4'^„ represents the numerical value of at point n of the grid cell shown in Fig. [2l The error term sym- 
bolizes possible terms of the form cx /I'^Alnft,, ft-^A^ln/i, /I'^AlnA and /i^A^lnA. In estimating this discretization 
error we must recall that the field is continuous, but not diffcrentiable — we expect derivatives of to diverge 
logarithmically near the worldline. The various error terms from the above discretization scheme are expected to 
be proportional to /I'^A (or /i^A^), times derivatives of somewhere in the cell. Since these derivatives diverge 
logarithmically for /i, A ^ 0, we must allow for the above logarithmic error terms. Note, finally, that we have dropped 
the contribution from the term cx coiO'^Yi.e in DIp^'r: This term is already at least of ©(A In A) near the particle 
[since cot ~ —{0 — 7r/2) near 9 = 7r/2], and so its volume integral can be absorbed in the error term in Eq. ([67| . 

We next turn to evaluate the RHS of Eq. ([66]) . We use the asymptotic expansion ([53| . recalling Z^^ = —{fr/AjS^. 
The expansion terms specified in Eq. (j66|) are sufficient for our purpose, since the volume integral of the ©(p* Inp*) er- 
ror term contributes only at order 0{h'^ In h) (and higher), which is the order of error allowed for on the LHS of Eq. (|66p . 

In terms of the local polar coordinates p*,4>* defined in Eq. ([50]). the volume element is dV^ = 2/Q~^^^rQ ^p*d/9* dtd0*, 
and the volume integral reads 

/ Z^dV^ ^e~*" / dp,dtd<l>,rf[a{cl>,) + Pinp,ln{pj4)+p,r{^*)]+0{hHnh), (68) 

J cell 47rro/o -^ceii 

where tc is the value of t at the central point c. Here we applied a "constant phase" approximation, which is valid since, 
for m not too large, the factor e~™"* varies negligibly across one grid cell. We next expand the factor rf = r — 2M 
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1 II 

in the integrand as r/ = rp/o + 5r = tq/o + /o^*^* = J^o/o + /o P* cos(/)*, with higher-order terms of the expansion 
absorbed in the error term 0(/i^lnft,) of the integral. We observe, recaUing a((/),) cx cos sin^ , that the term 
cx (7'o/o)q;(0*) vanishes upon integrating d(/)*. All 0*-dependent terms in {rofo)P"^{4>*) similarly yield a vanishing 
contribution. Omitting higher-order terms, we are left with 

ZndV = -^e-™""*^ / dp, dtd(f>,p, \a{,p,) cos,/), + ro/o^' (An ln(/5,/4) + /Jq™)! . (69) 



coll 47rro J coll 

This integral can be evaluated explicitly — most easily by first transforming back to local Cartesian coordinates Sr^, 60. 
We obtain 



cell 



ZndV = ^^ro/oe-™"^*^ {ao + ll/3in - 6/3o™ - (ao + 2An) [3C arctan(l/C) + (1/C) arctanC] 
z47r 



where 



+ 2 [2C2ao + A„(C' - 3)] In (cVl + C^) - 6A„ln (roA/(8P^/2)^ | 
^ 2A/.2zS',(t,), (70) 

8(1 - M/rp) f-i/2^A/M ^71^ 

Finally, equating Eqs. (p7|) and ([70]) and solving for ^'J^^, we obtain the finite-difference formula for worldline points: 

^r™ = [RHS of Eq. (Ell), with *™ *'^„] + /i^^g'^ (on the worldhne), (72) 

with a local error of 0{h^ In /i). 

We can summarize our finite-difference scheme for any grid point as follows. Given the values of the numerical field 
at points 2-8 of the grid cell depicted in Fig. [21 and assuming these values have already been adjusted as either all 
'in' or all 'out' (as described above), then the value of the field at point 1 is approximated by 

( +0{h*), point 1 outside T, 

= [RHS of Eq. dm)] -I- < /t^^Rc +0{h^), point 1 inside T, off worldline, (73) 
[ h'^Z^c +0{h^lnh), point 1 inside T, on worldline. 

The global (accumulated) error from points off the worldline is expected to be 0{h?). The error from worldline points 
accumulates only along the worldline (the number of points contributing to it scales as ~ ^/h), and is expected to 
dominate the global evolution error, with contribution of 0{h? In/i). We thus expect our scheme to converge at least 
linearly, but note that the above logarithmic error terms are likely to deter quadratic convergence. 

How could one eliminate the dominant logarithmic error terms, in order to assure quadratic convergence? The 
occurrence of such terms can be traced back to the logarithmic divergence of the R-field derivatives at the worldline, 
which is a feature of the leading-order puncture scheme adopted here. A natural solution to the problem could 
be offered within a higher-order puncture scheme, which incorporates a differentiable (C^) R-field. We leave the 
formulation of such a scheme for future work. 



D. Boundary and initial conditions 

At the poles {9 — 0,7r) we apply boundary conditions as in the vacuum case — see Eq. ([551) . initial conditions, 
we simply set the field (both ^E" and 5'r) to zero along u = uq and v = vq. This produces a burst of spurious radiation, 
mainly at the particle's initial location and at the intersection of dT with the initial surfaces. The spurious waves die 
off at late time, gradually unveiling the physically-meaningful field. In application of the code, one must monitor the 
effect of residual spurious waves, by testing the stationarity of the late-time field. 



VI. RESULTS AND CODE VALIDATION 



A. Sample results 



We present results for several sample cases, highlighting a few generic features of the numerical solutions generated 
using the above puncture scheme. The plots in Fig. [3 show numerical results for the modes m = 0, 1, 2 of the scalar 
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FIG. 7: Sample numerical results for ro = 7M (r*o = 8.8326M), showing each of the 3 modes ^'"-"'i-^ along 3 different 
slice-cuts of the 2+1-D domain: The upper left figure shows |4'™(r,)| a,t 9 = ■k/2 and t = 500M; the upper right figure 
shows |*'"(6')| at r = 7M and t = 500M; and the lower figure shows |*R(t)| at r = 7M and 9 = tt/2 (namely, along the 
particle's worldline). The dimensions of the auxiliary worldtube here are Sr, = 7.5M and Se — n/A. In the two spatial slices 
we display the full field outside the worldtube, and the residual field = - *p inside it. (Recall *p is the puncture 
function, given analytically.) Inside the worldtube we also indicate, in dotted line, the full (divergent) field, obtained through 
_ _|_ ^j^^ Spatial sUccs, dots along the graphs mark the location of actual numerical grid points. The various 

features of these solutions are discussed in the text. 



field for a circular geodesic orbit with radius ro = 7M {r^,o — 8.8326M). All solutions were obtained with a 
grid resolution oi h = M/4 and A = n/40, giving a ratio A/h ^ 0.31M^^ just above the Courant limit. (Unlike 
in the vacuum case, near the particle the resolution requirement in the longitudunal direction is as high as in the 
radial direction, which requires us to lower the ratio A/h.) The worldtube dimensions were taken as 5r, = 7.5M and 
6$ = 7r/4. The evolution starts at t = and ends at t = lOOOM (covering roughly 8 orbital periods). 

We highlight a few of the features visible in these plots: (i) The early stage of the numerical evolution is dominated 
by noise from spurious initial waves. These die off within ^ 1 orbital period, giving way to stationary behavior at late 
time, (ii) The full field calculated outside the worldtube merges smoothly with + \I'p across the boundaries of 
the worldline, as expected, (iii) The numerical variable is continuous at the particle, and has a well-defined value 
there. This, of course, makes it miich more tractable numerically than the original, divergent field 'J'™, (iv) The value 
of ^'p' at the particle drops rapidly with increasing m. (v) The residual field is asymmetric about the particle in 
the radial direction, reflecting the slight anisotropy of the curved background spacetime. It is this asymmetry in the 
field that gives rise to the SF effect. A scheme for constructing the physical SF from will be presented elsewhere. 
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B. Tests of code 

In what follows we demonstrate the numerical robustness of our code by (i) demonstrating point-wise numerical 
convergence, and (ii) showing that our numerical solutions depend only weakly on the dimensions of the auxiliary 
worldtube, and that this dependence gets ever weaker with improving resolution. We then validate our numerical 
solutions (and the entire puncture scheme) by comparing with results obtained from our 1+1-D companion code. 

1. Numerical convergence 

As in the vacuum case, we examined the point-wise convergence of our scheme by comparing solutions obtained 
with different grid resolutions. Figure [5] demonstrates the convergence properties of our solutions, for tq — 6.1M and 
m = 0, 1. The following features are manifest: (i) Away from the particle, and after the decay of the initial spurious 
waves, the numerical solutions show an approximate quadratic convergence, (ii) Near the particle (within a distance 
of a few M) the convergence is not uniform, and more difficult to characterize — but appears to be better than linear 
everywhere, (iii) On the particle itself, the convergence seems, once again, slower than quadratic and faster than 
linear. We have confirmed the generality of these features by experimenting with a range of orbital radii ro and modes 
m, and different auxiliary worldtube dimensions. We suspect that what slows down the convergence near the particle 
are the logarithmic error terms discussed above. Possible ways to improve the convergence of the scheme in future 
work will be discussed below. 



2. Dependence on worldtube dimensions 

It is important to establish that our numerical solutions do not depend on the dimensions of the auxiliary wolrdtube 
(modulo discretization error). We have tested the code with various worldtube dimensions, and present typical results 
in Figs.[9land llOI The figures compare between solutions obtained using two different worldtubes: One with dimensions 
{6r, , Sg) = (1.25A/, 7r/4), and the other with dimensions {Sr, , Sg) = (2.5Af, 7r/2). Evidently, the value of the calculated 
field is only very slightly affected by the different choice of worldtube, and the tiny differences appear to diminish 
rapidly with improving resolution. Figs. and 1101 demonstrate this behavior for ro ~ 7M and m = 0, 1, but we 
observe similar behavior for other ro and m. 

3. Comparison with 1+1-D solutions 

A good quantitative test of our code is provided by comparing the 2-1- 1-D solutions with solutions obtained using 
a H-l-D evolution code. To obtain 1-l-l-D solutions for a scalar charge in a circular orbits, we extended our vacuum 
1-l-l-D code to incorporate a source particle, using the prescription of Ref. [3l[. The code calculates individual 
multipole modes m of the scalar field. To allow comparison with the 2-l-l-D code, we decomposed the 2-l-l-D 
numerical solutions into their individual € modes (by integrating numerically with respect to against individual 
Legendre functions with given m). Results from such comparison, for rg = IM and m = 0, 1, 2, are shown in Fig. 
rm In all cases examined we find convincing agreement between the 2-l-l-D and 1-l-l-D solutions. 

VII. SUMMARY AND DISCUSSION 

In this work we began developing the computational framework to facilitate evolution of black hole perturbations 
from point particles in 2-l-l-D. This is mainly motivated within the context of SF calculations in Kerr: Knowledge of 
the perturbation field near the particle is an essential input for any calculation of the local SF. The local field near 
a particle in Kerr orbits has been studied so far mainly in a 1-l-l-D framework (i.e., through a spherical-harmonic 
decomposition), and so we started our analysis by exploring the field behavior in 2-l-l-D. Specializing to a scalar 
field, we established that each azimuthal rn-mode of the perturbation generically shows a logarithmic divergence 
near the particle. We then devised a numerical evolution scheme for the scalar field, based on approximating the 
divergent piece of the field analytically, and solving for the finite (and continuous) residual part. We demonstrated 
the applicability of this "puncture" scheme in the test case of circular orbits in Schwarzschild (but working in 2-l-l-D). 
For this, we developed a new 2-l-l-D evolution code, which we tested for numerical robustness and by demonstrating 
agreement with solutions obtained using other methods. We found that the scheme successfully resolves the residual. 
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FIG. 8: Numerical convergence test for the puncture scheme. Here we set tq = 6.1M (r»o — 7.5357M), and take worldtube 
dimensions of Sr, = 5M, 5e = n/3. We plot the relative differences (5*;;\ = \ (*^) - - | , forn = 2, 4, where 

^'(^) is the solution obtained with resolution h = (10M/7r)A = M/n. Outside the worldtube we show the relative differences in 
the full field and inside the worldtube— in the regular variable *r. Upper-left panel: (5*™i(r.) for {1,6) = (f OOOAf, 7r/2) 
(crossing the particle). Upper-right panel: 5'^^i{9) for = (lOOOM, ro) (crossing the particle), -left panel: 5'i!^i(9) 

for (t,r) = (I000M,4.2766A/) (off the particle). Lower-nght panel: (5*™i(t) for {e,r) = (7r/2,ro), i.e., along the particle's 
worldline. Away from the particle, the convergence is approximately quadratic; near the particle the convergence rate is harder 
to characterize, but remains at least linear. 



sub-dominant behavior of the scalar field near the particle. Whether our code, in its current form (and with realistic 
numerical resolution), allows sufficient accuracy for precise SF calculations remains to be explored, and we leave this 
for future work. 

The main strength of our time-domain approach is in the fact that it is now rather straightforward to extend the 
analysis to problems which are more astrophysically interesting. Our scalar field code can be readily extended to deal 
with inclined and eccentric orbits, and generalization to Kerr spacetime could be achieved rather straightforwardly 
based on the existing platform. We envisage applying the same numerical approach for solving the gravitational 
perturbation equations (in the Lorenz gauge), but we anticipate this would require much preparatory formulation 
work (to cast the equations in a form suitable for 2-f 1-D evolution). 

Before attempting further extensions/applications of the code, one may consider a few possible improvements of 
the numerical method. Firstly, working with a higher-order puncture scheme could prove very beneficial. In our 
leading-order scheme, the regularized field is not differentiable at the location of the particle (derivatives of ^'^ 
diverge there logarithmically, in a direction-dependent manner), which complicates the analysis. This can be cured 
by including higher-order terms in the definition of the analytic puncture \E'p. A second-order puncture can readily 
be constructed based on Eq. ((5^ above, which should yield a differentiable field (with nearly no extra cost 
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FIG. 9: Independence of the numerical solutions on the choice of worldtube dimensions — demonstrated here for ro = 7M, 
m = 0, 1. The graphs compare the solutions obtained using two different auxiliary worldtubes: one with dimensions , 5e) = 
(1.25M, 7r/4) (dashed line), and the other with dimensions {5r, , Se) ~ (2.5A/, n/2) (solid line). The dotted line is the full solution 
+^p', as obtained with the larger worldtube. The left panel displays the behavior as a function of r* at {t, 9) = (500M, 7r/2), 
and the right panel shows the behavior as a function of 9 at (t,r) = (500M, ro). The two calculations agree well on the value 
of ^B. inside the small worldtube, and on the value of elsewhere. It is demonstrated in Fig. [TO] below (for the m = case) 
that the tiny discrepancy between the two solutions tends to zero with increasing grid resolution. 
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FIG. 10: Data corresponding to the two upper plots of Fig. [9] Shown here, for m = 0, is the ratio between the two numerical 
solutions obtained with different worldtubes (i.e., the ratio between the dash and solid lines in Fig.[2}. The ratio is displayed for 
three different grid resolutions [h = M/4, M/8, M/16, where in each case A = {n/10)M-^h]. The smaU discrepancy between 
the two solutions diminishes rapidly with increasing resolution, suggesting that our numerical solutions are insensitive to the 
choice of auxiliary worldtube at the continuum limit, as should be expected. 



in computation time). This should have the following benefits: (i) The source term in the i?-field equation would 
no longer have a strong, divergence near the particle — it would instead diverge logarithmically, which is much 
easier to accommodate numerically, (ii) No logarithmic error terms of the sort discussed in Sec. IV CI would occur, 
which should settle the quadratic convergence of the scheme near the particle, (iii) SF calculations require the field's 
derivatives at the particle. The i?-field calculated using our leading-order puncture would therefore require further 
regularization, whereas, at least in principle, the SF should be accessible directly from a differentiable i?-filed coming 
from a second-order puncture scheme (cf. discussion below). 

Another possible improvement concerns the choice of initial data for the evolution. Expediting the damping of 
the spurious initial waves would allow shorter runs and save in computational cost. An improved scheme could 
incorporate smooth approximate initial data. Another idea, easier to implement, is to use interpolated solutions 
from low-resolution runs as initial conditions for high-resolution evolution. We are planning to incorporate the latter 
scheme in future applications of our code. 

Finally, we briefly discuss the application of our method to SF calculations. The standard "mode-sum" formula 
for the SF in Kerr [7| requires as input the individual £, m multipole modes of the perturbation field. To apply the 
mode-sum formula in its standard form, we would therefore have to decompose our numerical m-mode solutions into 
their individual ^-mode components. A more direct approach would be to access the SF directly from our regular 
m-mode fields \E'p^. We are currently formulating such an "m-mode sum" scheme for the SF, and will present it 



22 




FIG. 11: Left panel: Comparison of 2+1-D solutions (dashed line) with 1+1-D solutions (solid line), for individual multipole 
modes {£,m) — (0,0), (1, 1), (2,2). The 2+1-D modes were obtained by decomposing the original 2+1-D solutions into their 
individual £ components (using numerical integration). The data here is for ro = 7M, and the field is extracted at t = lOOOM 
in both codes. Right panel: The relative difference between the 2+1-D and 1+1-D solutions shown on the left. All modes 
compare to within less than 1% difference (or far better), providing a strong validation test for the 2+1-D code. 

elsewhere [s^. The proposed scheme could use either the first-order puncture solutions calculated here (which 
would then require us to apply a certain local-averaging procedure), or, in its simpler form, it could use the solutions 
of a future second-order puncture code. 
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APPENDIX A: TABLES OF POLYNOMIALS 

We tabulate here the various polynomials appearing in the expressions for $™ [Eq. and [Eq. iP5|) with 

Eq. pg)) ]. for the 6 modes m ~ 0-5. The polynomials and are listed in Tablc[Tl The polynomials and p™^ 
(with n — 1-4) are listed in Table iHl 

APPENDIX B: THE LOCAL SOURCE COEFFICIENTS f3^ 

The coefficients /3™ in Eq. ([55)1 are given by 

(ro - 2My^P^^ 

where are dimensionless polynomials in M /tq. A list of these polynomials, for m = 0-5, is provided in Table Hill 
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TABLE I: The Polynomials and appearing in Eq. (|42p . These are given in terms of the dimensionless quantity p = 
p/(2P^/^). 
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A procedure to obtain the SF directly from the Weyl scalars through direct regularization of the latter has been suggested 
recently, but has not yet been implemented for orbiting particles, 
[34] One may consider an alternative puncture function, obtained by replacing 2(1 — cos Sip) sin^ 5(p in Eq. p8|l . This has the 
advantage that the odd-m modes of "I>p and Sr can be written in terms of elementary functions. The disadvantage is that, 
to avoid the singularity of this alternative $p a.t p = 0, S(p — n, one has to introduce a cut-off on $p at some \S(p\ < tt, 
which then generates distributional contributions to the source modes Sr, complicating their form considerably. 
[35] This can be seen by considering the volume integral of □4'r over a small 3-ball (in E) surrounding the particle, at the 
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TABLE III: The dimensionless polynomials /3™ appearing in Eq. ()B1|) . In this table C, = M /tq. 



limit where the radius of the ball tends to zero. Using the Gauss theorem, this can be converted to a surface integral of 
'I'R.ci over the 2-sphere. By virtue of Eq. (|35|l we have that "I>r is bounded at the particle, and that the gradient ^r.q can 
at most diverge as ~ 1/eo there. Hence, the surface integral of 'I>R,a vanishes as the radius of the 2-sphere is taken to zero, 
implying □$r (and hence also S'r) contains no Dirac deltas. 
] The relation between local and global errors in the scheme (|61|l can be explained as follows: At each grid point, the field 
accumulates local errors from oc h"^ points belonging to the same S=const slice. [Note that the leading-order contribution 
to in Eq. (|6ip comes from points 2-4, which lie on the same &=const slice as point 1.] Assuming the local 0{h'^) errors 
are not strongly correlated, they accumulate to give a global error 



